#1) Study 1: Poland and Spain ####

#1.0) Load data ####

ples <- dplyr::left_join(readRDS("./data/exp1.RDS"),
                         readRDS("./data/respondent.RDS") %>% 
                           select(id, lrecon, income_cat, income, NatInt),
                         by = "id") %>% 
  mutate(individual = individual / 100,
         societal = societal / 100,
         foreign = foreign / 100,
         isdistance = individual - societal,
         ifdistance = individual - foreign,
         sfdistance = societal - foreign,
         is_dist = abs(isdistance),
         if_dist = abs(ifdistance),
         sf_dist = abs(sfdistance),
         is_dist_bi = factor(ifelse(isdistance > 0, "Advantageous", "Disadvantageous"),
                             levels = c("Advantageous", "Disadvantageous")),
         if_dist_bi = factor(ifelse(ifdistance > 0, "Advantageous", "Disadvantageous"),
                             levels = c("Advantageous", "Disadvantageous")),
         sf_dist_bi = factor(ifelse(sfdistance > 0, "Advantageous", "Disadvantageous"),
                             levels = c("Advantageous", "Disadvantageous")),
         iweight = ifelse(individual < 0, (1/(1-(10/41)^2)), 1),
         sweight = ifelse(societal < 0, (1/(1-(10/41)^2)), 1),
         fweight = ifelse(foreign < 0, (1/(1-(10/41)^2)), 1),
         isweight = ifelse(individual < 0 & societal < 0, (1/(1-(10/41))), 1),
         ifweight = ifelse(individual < 0 & foreign < 0, (1/(1-(10/41))), 1),
         sfweight = ifelse(foreign < 0 & societal < 0, (1/(1-(10/41))), 1),
         weight = ifelse(individual < 0 & societal >= 0 & foreign >= 0, iweight,
                         ifelse(individual < 0 & societal < 0 & foreign >= 0, isweight,
                                ifelse(individual < 0 & societal >= 0 & foreign < 0, ifweight,
                                       ifelse(individual >= 0 & societal < 0 & foreign >= 0, sweight,
                                              ifelse(individual >= 0 & societal < 0 & foreign < 0, sfweight,
                                                     ifelse(individual >= 0 & societal >= 0 & foreign < 0, fweight, 1)))))))


tapply(ples$rating, ples$country, mean, na.rm = T)
prop.table(table(as.numeric(haven::as_factor(ples$rating))>4, ples$country),2)

weighted.mean(ples$individual, ples$weight)
weighted.mean(ples$societal, ples$weight)
weighted.mean(ples$foreign, ples$weight)

ples <- ples %>% 
  mutate(is_dist_cat = cut(is_dist, quantile(is_dist),
                           labels = c("Q1", "Q2", "Q3", "Q4")),
         is_dist_cat = ifelse(is.na(is_dist_cat), "1", is_dist_cat),
         if_dist_cat = cut(if_dist, quantile(if_dist),
                           labels = c("Q1", "Q2", "Q3", "Q4")),
         if_dist_cat = ifelse(is.na(if_dist_cat), "1", if_dist_cat),
         sf_dist_cat = cut(sf_dist, quantile(sf_dist),
                           labels = c("Q1", "Q2", "Q3", "Q4")),
         sf_dist_cat = ifelse(is.na(sf_dist_cat), "1", sf_dist_cat)
  )
table(ples$is_dist_cat)
table(ples$if_dist_cat)
table(ples$sf_dist_cat)

ples$is_dist_cat <- factor(ples$is_dist_cat,
                           labels = c("Q1", "Q2", "Q3", "Q4"))

ples$if_dist_cat <- factor(ples$if_dist_cat,
                           labels = c("Q1", "Q2", "Q3", "Q4"))

ples$sf_dist_cat <- factor(ples$sf_dist_cat,
                           labels = c("Q1", "Q2", "Q3", "Q4"))

#1.1) Model 0: Descriptive analysis ####

plesm0r <- lm(rating ~ individual + societal + foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
lmtest::coeftest(plesm0r, vcov = sandwich::vcovCL, cluster = ~id)

plesm0c <- lm(choice ~ individual + societal + foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
lmtest::coeftest(plesm0c, vcov = sandwich::vcovCL, cluster = ~id)

texreg::screenreg(list(prepTex(plesm0r, data = ples), prepTex(plesm0c, data = ples)))

marginaleffects::avg_comparisons(plesm0r, variable = list(individual=c(0, 4), societal=c(0, 4), foreign=c(0, 4)), vcov = sandwich::vcovCL) #Attention: ranges used different across dimensions, see paper
marginaleffects::avg_comparisons(plesm0c, variable = list(individual=c(0, 4), societal=c(0, 4), foreign=c(0, 4)), vcov = sandwich::vcovCL) #Attention: ranges used different across dimensions, see paper


#1.1.1) Table 2 ####

texreg::texreg(list(prepTex(plesm0r, data = ples), prepTex(plesm0c, data = ples)),
               file = "./tables/Table2.tex",
               custom.model.names = c("Rating", "Choice"),
               custom.coef.map = list("individual" = "Individual gains",
                                      "societal" = "Societal gains",
                                      "foreign" = "Foreign gains",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Year of Implementation (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               caption = "Individual, societal, and foreign gains and PTA support",
               caption.above = T,
               label = "tab:table2", 
               custom.note = "\\parbox{0.6\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on respondents. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities. The gains/losses were divided by 100 to make the coefficients easier to interpret. They thus range from -1 to + 3.}",
               use.packages = F)


#1.2) Model 1: Investigation of Absolute Distance ####

plesm1sr <- lm(rating ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sc <- lm(choice ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fr <- lm(rating ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fc <- lm(choice ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)

texreg::screenreg(list(prepTex(plesm1sr, data = ples), prepTex(plesm1sc, data = ples),
                       prepTex(plesm1fr, data = ples), prepTex(plesm1fc, data = ples)))

marginaleffects::avg_comparisons(plesm1sr, variable  = list(is_dist=c(0, 4)), vcov = sandwich::vcovCL)
marginaleffects::avg_comparisons(plesm1sc, variable = list(is_dist=c(0, 4)), vcov = sandwich::vcovCL)
marginaleffects::avg_comparisons(plesm1fr, variable = list(if_dist=c(0, 4)), vcov = sandwich::vcovCL)
marginaleffects::avg_comparisons(plesm1fc, variable = list(if_dist=c(0, 4)), vcov = sandwich::vcovCL)

plesm1sr_cat <- lm(rating ~ is_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sc_cat <- lm(choice ~ is_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fr_cat <- lm(rating ~ if_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fc_cat <- lm(choice ~ if_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)


#1.2.1) Table D1 ####
texreg::texreg(list(prepTex(plesm1sr, data = ples), prepTex(plesm1sc, data = ples),
                    prepTex(plesm1fr, data = ples), prepTex(plesm1fc, data = ples)),
               file = "./tables/TableD1.tex",
               custom.model.names = c("$I-S$ (Rating)", "$I-S$ (Choice)", "$I-F$ (Rating)", "$I-F$ (Choice)"),
               custom.coef.map = list("is_dist" = "Abs. distance",
                                      "if_dist" = "Abs. distance",
                                      "individual" = "Individual gains",
                                      "societal" = "Societal gains",
                                      "foreign" = "Foreign gains",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = .85,
               caption = "Absolute distance and PTA support",
               caption.above = T,
               label = "tab:tableD1", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on respondents. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.3 Model 2: Absolute Distance and (Dis-)Advantageous Inequality ####

plesm2sr <- lm(rating ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm2sc <- lm(choice ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm2fr <- lm(rating ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm2fc <- lm(choice ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)

texreg::screenreg(list(prepTex(plesm2sr, data = ples), prepTex(plesm2sc, data = ples),
                       prepTex(plesm2fr, data = ples), prepTex(plesm2fc, data = ples)))

marginaleffects::avg_comparisons(plesm2sr, variable  = list(is_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::avg_comparisons(plesm2sc, variable = list(is_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::avg_comparisons(plesm2fr, variable = list(if_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "if_dist_bi")
marginaleffects::avg_comparisons(plesm2fc, variable = list(if_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "if_dist_bi")


#1.3.1) Table D2 ####

texreg::texreg(list(prepTex(plesm2sr, data = ples), prepTex(plesm2sc, data = ples),
                    prepTex(plesm2fr, data = ples), prepTex(plesm2fc, data = ples)),
               file = "./tables/TableD2.tex",
               custom.model.names = c("$I-S$ (Rating)", "$I-S$ (Choice)", "$I-F$ (Rating)", "$I-F$ (Choice)"),
               custom.coef.map = list("is_dist" = "Abs. distance",
                                      "if_dist" = "Abs. distance",
                                      "is_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "if_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "is_dist:is_dist_biDisadvantageous" = "Abs. dist $\\times$ Disadv. ineq.",
                                      "if_dist:if_dist_biDisadvantageous" = "Abs. dist $\\times$ Disadv. ineq.",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = 0.85,
               caption = "Absolute distance, (dis-)advantageous inequality and PTA support",
               caption.above = T,
               label = "tab:tableD2", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on individuals. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.4 Model 3: Categorical absolute Distance and (Dis-)Advantageous Inequality ####


plesm3sr <- lm(rating ~ is_dist_cat * is_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm3sc <- lm(choice ~ is_dist_cat * is_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm3fr <- lm(rating ~ if_dist_cat * if_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm3fc <- lm(choice ~ if_dist_cat * if_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)

texreg::screenreg(list(prepTex(plesm3sr, data = ples), prepTex(plesm3sc, data = ples),
                       prepTex(plesm3fr, data = ples), prepTex(plesm3fc, data = ples)))

marginaleffects::slopes(plesm3sr, variable = "is_dist_cat", vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::slopes(plesm3sc, variable = "is_dist_cat", vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::slopes(plesm3fr, variable = "if_dist_cat", vcov = sandwich::vcovCL, by = "if_dist_bi")
marginaleffects::slopes(plesm3fc, variable = "if_dist_cat", vcov = sandwich::vcovCL, by = "if_dist_bi")


#1.4.1) Table D4 ####

texreg::texreg(list(prepTex(plesm3sr, data = ples), prepTex(plesm3sc, data = ples),
                    prepTex(plesm3fr, data = ples), prepTex(plesm3fc, data = ples)),
               file = "./tables/TableD4.tex",
               custom.model.names = c("$I-S$ (Rating)", "$I-S$ (Choice)", "$I-F$ (Rating)", "$I-F$ (Choice)"),
               custom.coef.map = list("is_dist_catQ2" = "Abs. dist. (Q2)",
                                      "is_dist_catQ3" = "Abs. dist. (Q3)",
                                      "is_dist_catQ4" = "Abs. dist. (Q4)",
                                      "if_dist_catQ2" = "Abs. dist. (Q2)",
                                      "if_dist_catQ3" = "Abs. dist. (Q3)",
                                      "if_dist_catQ4" = "Abs. dist.(Q4)",
                                      "is_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "if_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "is_dist_catQ2:is_dist_biDisadvantageous" = "Abs. dist. (Q2) $\\times$ Disadv. ineq.",
                                      "is_dist_catQ3:is_dist_biDisadvantageous" = "Abs. dist. (Q3) $\\times$ Disadv. ineq.",
                                      "is_dist_catQ4:is_dist_biDisadvantageous" = "Abs. dist. (Q4) $\\times$ Disadv. ineq.",
                                      "if_dist_catQ2:if_dist_biDisadvantageous" = "Abs. dist. (Q2) $\\times$ Disadv. ineq.",
                                      "if_dist_catQ3:if_dist_biDisadvantageous" = "Abs. dist. (Q3) $\\times$ Disadv. ineq.",
                                      "if_dist_catQ4:if_dist_biDisadvantageous" = "Abs. dist. (Q4) $\\times$ Disadv. ineq.",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = 0.85,
               caption = "Categorized absolute distance, (dis-)advantageous inequality and PTA support",
               caption.above = T,
               label = "tab:tableD4", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on individuals. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)


#1.5) Preparation Figures 1 and E1 ####

#Plot1
predm1sr <- marginaleffects::predictions(plesm1sr, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1)),
                                         vcov = sandwich::vcovCL(plesm1sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist) %>% 
  mutate(xvar = "IS",
         dv = "R")

predm1sc <- marginaleffects::predictions(plesm1sc, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1)),
                                         vcov = sandwich::vcovCL(plesm1sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist) %>% 
  mutate(xvar = "IS",
         dv = "C")

predm1fr <- marginaleffects::predictions(plesm1fr, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1)),
                                         vcov = sandwich::vcovCL(plesm1fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist) %>% 
  mutate(xvar = "IF",
         dv = "R")

predm1fc <- marginaleffects::predictions(plesm1fc, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1)),
                                         vcov = sandwich::vcovCL(plesm1fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist) %>% 
  mutate(xvar = "IF",
         dv = "C")


predm1sr_cat <- marginaleffects::predictions(plesm1sr_cat, newdata = datagrid(is_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                             vcov = sandwich::vcovCL(plesm1sr_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist_cat) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "1")

predm1sc_cat <- marginaleffects::predictions(plesm1sc_cat, newdata = datagrid(is_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                             vcov = sandwich::vcovCL(plesm1sc_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist_cat) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "1")

predm1fr_cat <- marginaleffects::predictions(plesm1fr_cat, newdata = datagrid(if_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                             vcov = sandwich::vcovCL(plesm1fr_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist_cat) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "1")

predm1fc_cat <- marginaleffects::predictions(plesm1fc_cat, newdata = datagrid(if_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                             vcov = sandwich::vcovCL(plesm1fc_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist_cat) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "1")

axis_pos <- rbind(ples %>% 
                    group_by(is_dist_cat) %>% 
                    summarise(pos = mean(is_dist)) %>% 
                    filter(!is.na(is_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = is_dist_cat) %>% 
                    mutate(xvar = "Individual - Societal"),
                  ples %>% 
                    group_by(if_dist_cat) %>% 
                    summarise(pos = mean(if_dist)) %>% 
                    filter(!is.na(if_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = if_dist_cat) %>% 
                    mutate(xvar = "Individual - Foreign")) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("Individual - Societal", "Individual - Foreign")))

bin_1 <- rbind(predm1sr_cat, predm1sc_cat, predm1fr_cat, predm1fc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1") %>% 
  filter(dv == "Rating") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar")) %>% 
  mutate(pos = pos * 100)


rbind(predm1sr, predm1sc, predm1fr, predm1fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  geom_pointrange(data = bin_1, aes(x = pos, y = estimate, ymin = lower, ymax=upper), size = 0.1) +
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p1

bin_1c <- rbind(predm1sr_cat, predm1sc_cat, predm1fr_cat, predm1fc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1") %>% 
  filter(dv == "Choice") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar")) %>% 
  mutate(pos = pos * 100)

rbind(predm1sr, predm1sc, predm1fr, predm1fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  geom_pointrange(data = bin_1c, aes(x = pos, y = estimate, ymin = lower, ymax=upper), size = 0.1) +
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support")  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p1c

#Plot 2
predm2sr <- marginaleffects::predictions(plesm2sr, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm2sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "2")

predm2sc <- marginaleffects::predictions(plesm2sc, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm2sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "2")

predm2fr <- marginaleffects::predictions(plesm2fr, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm2fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "2")

predm2fc <- marginaleffects::predictions(plesm2fc, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm2fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "2")


#bins for Plot 2#


predm3sr <- marginaleffects::predictions(plesm3sr, newdata = datagrid(is_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm3sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist_cat, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "3")

predm3sc <- marginaleffects::predictions(plesm3sc, newdata = datagrid(is_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm3sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist_cat, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "3")

predm3fr <- marginaleffects::predictions(plesm3fr, newdata = datagrid(if_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm3fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist_cat, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "3")

predm3fc <- marginaleffects::predictions(plesm3fc, newdata = datagrid(if_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                         vcov = sandwich::vcovCL(plesm3fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist_cat, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "3")

axis_pos <- rbind(ples %>% 
                    group_by(is_dist_cat) %>% 
                    summarise(pos = mean(is_dist)) %>% 
                    filter(!is.na(is_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = is_dist_cat) %>% 
                    mutate(xvar = "Individual - Societal"),
                  ples %>% 
                    group_by(if_dist_cat) %>% 
                    summarise(pos = mean(if_dist)) %>% 
                    filter(!is.na(if_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = if_dist_cat) %>% 
                    mutate(xvar = "Individual - Foreign")) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("Individual - Societal", "Individual - Foreign")),
         pos = pos * 100)

bin_3 <- rbind(predm3sr, predm3sc, predm3fr, predm3fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 3") %>% 
  filter(dv == "Rating") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar"))

rbind(predm2sr, predm2sc, predm2fr, predm2fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_pointrange(data = bin_3, aes(x = pos, y = estimate, ymin = lower, ymax=upper, colour = mlevel), size = 0.1, position = position_dodge(width = 10)) +
  geom_line() + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = c(0.4, 0.1))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p2

bin_3c <- rbind(predm3sr, predm3sc, predm3fr, predm3fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 3") %>% 
  filter(dv == "Choice") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar"))

rbind(predm2sr, predm2sc, predm2fr, predm2fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 2",
         xlevel = 100 * xlevel) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_pointrange(data = bin_3c, aes(x = pos, y = estimate, ymin = lower, ymax=upper, colour = mlevel), size = 0.1, position = position_dodge(width = 10)) +
  geom_line() + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = c(0.4, 0.1))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p2c

ggpubr::ggarrange(p1, p2)

ggsave(filename = "./figures/Figure1.pdf",
       width = 20, height = 9, units = "cm")

ggpubr::ggarrange(p1c, p2c)

ggsave(filename = "./figures/FigureE1.pdf",
       width = 20, height = 9, units = "cm")

# 1.6) Loss aversion ####

ples$ind_cat1 <- factor(cut(ples$individual, breaks=c(-1.1, -0.1, .9, 1.9, 3),
                            labels=c('neg.', 'small', 'medium', 'high'), levels=c('neg.', 'small', 'medium', 'high')))
ples$soc_cat1 <- factor(cut(ples$societal, breaks=c(-1.1, -0.1, .9, 1.9, 3),
                            labels=c('neg.', 'small', 'medium', 'high'), levels=c('neg.', 'small', 'medium', 'high')))
ples$for_cat1 <- factor(cut(ples$foreign, breaks=c(-1.1, -0.1, .9, 1.9, 3),
                            labels=c('neg.', 'small', 'medium', 'high'), levels=c('neg.', 'small', 'medium', 'high')))
table(ples$ind_cat1, ples$soc_cat1)

plesm4sr <- lm(rating ~ ind_cat1 * soc_cat1 + ind_cat1 * for_cat1 + trade_volume + partner_size + implementation + country, ples, weight = weight)

# 1.6.1) Figure E2 ####

bind_rows(avg_comparisons(plesm4sr, variables = list(soc_cat1=c("neg.", "high")), by="ind_cat1",
                          vcov = sandwich::vcovCL(plesm4sr, cluster = ~id)),
          avg_comparisons(plesm4sr, variables = list(for_cat1=c("neg.", "high")), by="ind_cat1",
                          vcov = sandwich::vcovCL(plesm4sr, cluster = ~id))) %>% 
  mutate(term = factor(ifelse(term == "soc_cat1", "Societal gains",
                              ifelse(term == "for_cat1", "Foreign gains", NA)),
                       levels = c("Societal gains", "Foreign gains")),
         ind_cat1 = forcats::fct_recode(ind_cat1, negative = "neg.")) %>% 
  ggplot(., aes(x = ind_cat1, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange() +
  facet_grid(~term) +
  labs(x = "Individual gains", y = "Marginal effect of high gains vs negative gains") +
  theme_minimal()  +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureE2.pdf",
       width = 20, height = 9, units = "cm")

#1.6.2) Models excluding losses

plesm1sra <- lm(rating ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples[ples$individual>=0 & ples$societal>=0, ], weight = weight)
plesm1sca <- lm(choice ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples[ples$individual>=0 & ples$societal>=0, ], weight = weight)
plesm1fra <- lm(rating ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples[ples$individual>=0 & ples$foreign>=0, ], weight = weight)
plesm1fca <- lm(choice ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples[ples$individual>=0 & ples$foreign>=0, ], weight = weight)

texreg::screenreg(list(prepTex(plesm1sra, data = ples[ples$individual>=0 & ples$societal>=0, ]), prepTex(plesm1sca, data = ples[ples$individual>=0 & ples$societal>=0, ]),
                       prepTex(plesm1fra, data = ples[ples$individual>=0 & ples$foreign>=0, ]), prepTex(plesm1fca, data = ples[ples$individual>=0 & ples$foreign>=0, ])))

# 1.6.2.1) Table D3 #### 

texreg::texreg(list(prepTex(plesm1sra, data = ples[ples$individual>=0 & ples$societal>=0, ]), prepTex(plesm1sca, data = ples[ples$individual>=0 & ples$societal>=0, ]),
                    prepTex(plesm1fra, data = ples[ples$individual>=0 & ples$foreign>=0, ]), prepTex(plesm1fca, data = ples[ples$individual>=0 & ples$foreign>=0, ])),
               file = "./tables/TableD3.tex",
               custom.model.names = c("$I-S$ (Rating)", "$I-S$ (Choice)", "$I-F$ (Rating)", "$I-F$ (Choice)"),
               custom.coef.map = list("is_dist" = "Abs. distance",
                                      "if_dist" = "Abs. distance",
                                      "individual" = "Individual gains",
                                      "societal" = "Societal gains",
                                      "foreign" = "Foreign gains",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = .85,
               caption = "Absolute distance and PTA support (no individual or societal/foreign losses)",
               caption.above = T,
               label = "tab:TableD3", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on respondents. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)


# 1.7) Results by country ####

plm2sr <- lm(rating ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation, ples, subset = country == "PL", weight = weight)
plm2sc <- lm(choice ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation, ples, subset = country == "PL", weight = weight)
plm2fr <- lm(rating ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation, ples, subset = country == "PL", weight = weight)
plm2fc <- lm(choice ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation, ples, subset = country == "PL", weight = weight)

esm2sr <- lm(rating ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation, ples,, subset = country == "ES", weight = weight)
esm2sc <- lm(choice ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation, ples,, subset = country == "ES", weight = weight)
esm2fr <- lm(rating ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation, ples,, subset = country == "ES", weight = weight)
esm2fc <- lm(choice ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation, ples,, subset = country == "ES", weight = weight)

# 1.7.1) Figure E3 ####

bind_rows(
  ples %>% 
    group_by(country) %>% 
    do(broom::tidy(lm(cbind(rating, choice) ~ is_dist * is_dist_bi + trade_volume + partner_size + implementation, data = ., weight = weight))) %>% 
    filter(grepl("is_dist", term)),
  ples %>% 
    group_by(country) %>% 
    do(broom::tidy(lm(cbind(rating, choice) ~ if_dist * if_dist_bi + trade_volume + partner_size + implementation, data = ., weight = weight))) %>% 
    filter(grepl("if_dist", term))) %>% 
  mutate(dim = factor(ifelse(grepl("is_", term), "Societal gains", "Foreign gains"),
                      levels = c("Societal gains", "Foreign gains")),
         term = factor(term,
                       levels = c("is_dist",
                                  "is_dist_biDisadvantageous",
                                  "is_dist:is_dist_biDisadvantageous",
                                  "if_dist",
                                  "if_dist_biDisadvantageous",
                                  "if_dist:if_dist_biDisadvantageous"),
                       labels = c("Abs. distance",
                                  "Disadv Ineq.",
                                  "Interaction",
                                  "Abs. distance",
                                  "Disadv Ineq.",
                                  "Interaction"))) %>% 
  ggplot(., aes(x =  term, y = estimate,
                ymin = estimate - 1.96 * std.error,
                ymax = estimate + 1.96 * std.error,
                shape = country,
                colour = country)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) +
  facet_grid(response ~ dim, scales = "free_y")  + 
  theme_minimal() + 
  labs(x = "Term", y = "Estimate") +
  scale_colour_discrete("Country") + 
  scale_shape_discrete("Country") + 
  theme(legend.position = c(0.6, 0.15)) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureE3.pdf",
       width = 20, height = 9, units = "cm")

#1.7.2) Figures E4 and E5 ####

plpredm2sr <- marginaleffects::predictions(plm2sr, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(plm2sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "2",
         country = "PL")

plpredm2sc <- marginaleffects::predictions(plm2sc, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(plm2sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "2",
         country = "PL")

plpredm2fr <- marginaleffects::predictions(plm2fr, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(plm2fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "2",
         country = "PL")

plpredm2fc <- marginaleffects::predictions(plm2fc, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(plm2fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "2",
         country = "PL")

espredm2sr <- marginaleffects::predictions(esm2sr, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(esm2sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "2",
         country = "ES")

espredm2sc <- marginaleffects::predictions(esm2sc, newdata = datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                      is_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(esm2sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "2",
         country = "ES")

espredm2fr <- marginaleffects::predictions(esm2fr, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(esm2fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "2",
         country = "ES")

espredm2fc <- marginaleffects::predictions(esm2fc, newdata = datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                      if_dist_bi = c("Disadvantageous", "Advantageous")),
                                           vcov = sandwich::vcovCL(esm2fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "2",
         country = "ES")

rbind(plpredm2sr, plpredm2fr, espredm2sr, espredm2fr) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         country = factor(country,
                          levels = c("PL", "ES"),
                          labels = c("Poland", "Spain")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 2",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  facet_grid(xvar~country, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = c(0.6, 0.15))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureE4.pdf",
       width = 20, height = 9, units = "cm")

rbind(plpredm2sc, plpredm2fc, espredm2sc, espredm2fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         country = factor(country,
                          levels = c("PL", "ES"),
                          labels = c("Poland", "Spain")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 2",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  facet_grid(xvar~country, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = c(0.6, 0.125))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureE5.pdf",
       width = 20, height = 9, units = "cm")

#1.8) Marginal effects for all figures ####

#1.8.1) Figures F1 and F2 ####

margm1sr <- marginaleffects::avg_slopes(plesm1sr, variable = "is_dist",
                                        vcov = sandwich::vcovCL(plesm1sr, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high) %>% 
  mutate(xvar = "IS",
         dv = "R")

margm1sc <- marginaleffects::avg_slopes(plesm1sc, variable = "is_dist",
                                        vcov = sandwich::vcovCL(plesm1sc, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high) %>% 
  mutate(xvar = "IS",
         dv = "C")

margm1fr <- marginaleffects::avg_slopes(plesm1fr, variable = "if_dist",
                                        vcov = sandwich::vcovCL(plesm1fr, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high) %>% 
  mutate(xvar = "IF",
         dv = "R")

margm1fc <- marginaleffects::avg_slopes(plesm1fc, variable = "if_dist",
                                        vcov = sandwich::vcovCL(plesm1fc, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high) %>% 
  mutate(xvar = "IF",
         dv = "C")


rbind(margm1sr, margm1sc, margm1fr, margm1fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(~dv, scales = "free_x", switch = "y") + 
  labs(x = "Absolute Distance", y = "Marg. effect") +
  theme(legend.position = "bottom") + 
  scale_colour_discrete("Comparison")  +
  scale_y_continuous(minor_breaks = NULL) -> p1m


#Plot 2
margm2sr <- marginaleffects::avg_slopes(plesm2sr, variable = "is_dist", by = "is_dist_bi",
                                        vcov = sandwich::vcovCL(plesm2sr, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "2")

margm2sc <- marginaleffects::avg_slopes(plesm2sc, variable = "is_dist", by = "is_dist_bi",
                                        vcov = sandwich::vcovCL(plesm2sc, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "2")

margm2fr <- marginaleffects::avg_slopes(plesm2fr, variable = "if_dist", by = "if_dist_bi",
                                        vcov = sandwich::vcovCL(plesm2fr, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "2")

margm2fc <- marginaleffects::avg_slopes(plesm2fc, variable = "if_dist", by = "if_dist_bi",
                                        vcov = sandwich::vcovCL(plesm2fc, cluster = ~id)) %>% 
  select(term, estimate, std.error, lower = conf.low, upper = conf.high, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "2")

rbind(margm2sr, margm2sc, margm2fr, margm2fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, group = mlevel, colour = mlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = c(0.3, 0.2)) + 
  scale_colour_discrete("Inequity") +
  scale_y_continuous(minor_breaks = NULL) -> p2m

ggpubr::ggarrange(p1m, p2m, widths = c(1.5,2))

ggsave(filename = "./figures/FigureF1.pdf",
       width = 20, height = 9, units = "cm")

rbind(margm1sr, margm1sc, margm1fr, margm1fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  #coord_flip() + 
  theme_minimal() +
  facet_grid(~dv, scales = "free_x", switch = "y") + 
  labs(x = "Absolute Distance", y = "Marg. effect") +
  theme(legend.position = "bottom") + 
  scale_colour_discrete("Comparison")   +
  scale_y_continuous(minor_breaks = NULL)-> p1mc

rbind(margm2sr, margm2sc, margm2fr, margm2fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, group = mlevel, colour = mlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = c(0.3, 0.2)) + 
  scale_colour_discrete("Inequity")   +
  scale_y_continuous(minor_breaks = NULL) -> p2mc

ggpubr::ggarrange(p1mc, p2mc, widths = c(1.5,2))

ggsave(filename = "./figures/FigureF2.pdf",
       width = 20, height = 9, units = "cm")

#1.8.2) Figures F3 and F4####

#plot1
margm1sr_cat <- marginaleffects::avg_slopes(plesm1sr_cat, variable = "is_dist_cat", 
                                            vcov = sandwich::vcovCL(plesm1sr_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "1",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm1sc_cat <- marginaleffects::avg_slopes(plesm1sc_cat, variable = "is_dist_cat",
                                            vcov = sandwich::vcovCL(plesm1sc_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "1",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm1fr_cat <- marginaleffects::avg_slopes(plesm1fr_cat, variable = "if_dist_cat", 
                                            vcov = sandwich::vcovCL(plesm1fr_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "1",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm1fc_cat <- marginaleffects::avg_slopes(plesm1fc_cat, variable = "if_dist_cat", 
                                            vcov = sandwich::vcovCL(plesm1fc_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "1",
         xlevel = substr(xlevel, start = 6, stop = 7))


rbind(margm1sr_cat, margm1sc_cat, margm1fr_cat, margm1fc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, colour = xlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = "none") + 
  scale_colour_discrete("Inequity")   +
  scale_y_continuous(minor_breaks = NULL) -> p3m


#plot2
margm3sr <- marginaleffects::slopes(plesm3sr, variable = "is_dist_cat", by = "is_dist_bi",
                                    vcov = sandwich::vcovCL(plesm3sr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "R",
         mod = "3",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm3sc <- marginaleffects::slopes(plesm3sc, variable = "is_dist_cat", by = "is_dist_bi",
                                    vcov = sandwich::vcovCL(plesm3sc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast, mlevel = is_dist_bi) %>% 
  mutate(xvar = "IS",
         dv = "C",
         mod = "3",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm3fr <- marginaleffects::slopes(plesm3fr, variable = "if_dist_cat", by = "if_dist_bi",
                                    vcov = sandwich::vcovCL(plesm3fr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "R",
         mod = "3",
         xlevel = substr(xlevel, start = 6, stop = 7))

margm3fc <- marginaleffects::slopes(plesm3fc, variable = "if_dist_cat", by = "if_dist_bi",
                                    vcov = sandwich::vcovCL(plesm3fc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = contrast, mlevel = if_dist_bi) %>% 
  mutate(xvar = "IF",
         dv = "C",
         mod = "3",
         xlevel = substr(xlevel, start = 6, stop = 7))

rbind(margm3sr, margm3sc, margm3fr, margm3fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance",
         mlevel = paste0(mlevel, "\ninequity")) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, colour = xlevel, shape = mlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = c(0.4, 0.2),
        legend.direction="horizontal") + 
  scale_colour_discrete("Group", guide = "none") + 
  scale_shape_discrete("") +
  annotate(geom="text", x=0.675, y=0.04, label="Q2",
           color=2) + 
  annotate(geom="text", x=1, y=0.02, label="Q3",
           color=3) +
  annotate(geom="text", x=1.325, y=-0.05, label="Q4",
           color=4)   +
  scale_y_continuous(minor_breaks = NULL) -> p4m

ggpubr::ggarrange(p3m, p4m, widths = c(1.5,2))

ggsave(filename = "./figures/FigureF3.pdf",
       width = 20, height = 9, units = "cm")


rbind(margm1sr_cat, margm1sc_cat, margm1fr_cat, margm1fc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance") %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, colour = xlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = "none") + 
  scale_colour_discrete("Inequity") +
  scale_y_continuous(minor_breaks = NULL) -> p3mc

rbind(margm3sr, margm3sc, margm3fr, margm3fc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         term = "Absolute Distance",
         mlevel = paste0(mlevel, "\ninequity")) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xvar, y = estimate, ymin = lower, ymax = upper, colour = xlevel, shape = mlevel)) + 
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) + 
  theme_minimal() +
  facet_grid(.~dv, scales = "free", switch = "y") + 
  labs(x = "", y = "") +
  theme(legend.position = c(0.4, 0.2),
        legend.direction="horizontal") + 
  scale_colour_discrete("Group", guide = "none") + 
  scale_shape_discrete("") +
  annotate(geom="text", x=0.7, y=0.03, label="Q2",
           color=2) + 
  annotate(geom="text", x=1, y=0.01, label="Q3",
           color=3) +
  annotate(geom="text", x=1.3, y=-0.05, label="Q4",
           color=4)   +
  scale_y_continuous(minor_breaks = NULL) -> p4mc

ggpubr::ggarrange(p3mc, p4mc, widths = c(1.5,2))

ggsave(filename = "./figures/FigureF4.pdf",
       width = 20, height = 9, units = "cm")

#1.9) Round effects ####

#1.9.1) Figure E6 ####

ples %>% 
  group_by(round) %>% 
  do(broom::tidy(lm(cbind(rating, choice) ~ individual + societal + foreign + trade_volume + partner_size + implementation + country, data = ., weight = weight))) %>% 
  filter(term %!in% c("(Intercept)", "countryPL")) %>% 
  mutate(term = factor(term,
                       levels = c("individual", "societal", "foreign",
                                  "trade_volumelarge", "partner_sizelarge economy", "implementation2027"),
                       labels = c("Individual", "Societal", "Foreign", 
                                  "Large trade\nvolume", "Large\npartner", "Implementation\nby 2027")),
         response = factor(response,
                           labels = c("Choice", "Rating"))) %>% 
  filter(term %in% c("Individual", "Societal", "Foreign")) %>% 
  ggplot(., aes(x = forcats::fct_rev(round), y = estimate, 
                ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) +
  coord_flip() +
  facet_grid(term~response, scales = "free_x") +
  labs(y = "Coefficient", x = "Round") + 
  theme_minimal()  +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureE6.pdf",
       width = 20, height = 13.5, units = "cm")


#1.9.2) Figure E7 ####

bind_rows(ples %>% 
            group_by(round) %>% 
            do(broom::tidy(lm(cbind(rating, choice) ~ is_dist + trade_volume + partner_size + implementation + country, data = ., weight = weight))) %>% 
            mutate(mod = "Societal"),
          ples %>% 
            group_by(round) %>% 
            do(broom::tidy(lm(cbind(rating, choice) ~ if_dist + trade_volume + partner_size + implementation + country, data = ., weight = weight))) %>% 
            mutate(mod = "Foreign")) %>% 
  filter(term %!in% c("(Intercept)", "countryPL")) %>% 
  mutate(term = factor(term,
                       levels = c("is_dist", "if_dist",
                                  "trade_volumelarge", "partner_sizelarge economy", "implementation2027"),
                       labels = c("Individual - Societal", "Individual - Foreign", 
                                  "Large trade\nvolume", "Large\npartner", "Implementation\nby 2027")),
         response = factor(response,
                           labels = c("Choice", "Rating"))) %>% 
  filter(term %in% c("Individual - Societal", "Individual - Foreign")) %>% 
  ggplot(., aes(x = forcats::fct_rev(round), y = estimate, 
                ymin = estimate - 1.96*std.error, ymax = estimate + 1.96*std.error)) +
  geom_hline(yintercept = 0, lty = "dotted") + 
  geom_pointrange(position = position_dodge(width = 1)) +
  coord_flip() +
  facet_grid(term~response, scales = "free_x") +
  theme_minimal() +
  labs(y = "Coefficient", x = "Round") + 
  scale_y_continuous(minor_breaks = NULL) +
  NULL

ggsave(filename = "./figures/FigureE7.pdf",
       width = 20, height = 9, units = "cm")

# 1.10) Remove respondents with inconsistent behaviour ####

ples <- dplyr::left_join(ples, 
                         dplyr::left_join(ples %>% 
                                            select(id, side, round, rating) %>% 
                                            pivot_wider(names_from = side, names_prefix = "rating_", values_from = rating),
                                          ples %>% 
                                            select(id, side, round, choice) %>% 
                                            pivot_wider(names_from = side, names_prefix = "choice_",values_from = choice),
                                          by = c("id", "round")) %>% 
                           mutate(keep = ifelse((rating_A >= rating_B) & (choice_A > choice_B), 1, 
                                                ifelse((rating_B >= rating_A) & (choice_B > choice_A), 1, 0))) %>% 
                           select(id, round, keep),
                         by = c("id", "round"))

prop.table(table(ples$keep))

plesm1srkeep <- lm(rating ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, subset = keep == 1, weight = weight)
plesm1sckeep <- lm(choice ~ is_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, subset = keep == 1, weight = weight)
plesm1frkeep <- lm(rating ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, subset = keep == 1, weight = weight)
plesm1fckeep <- lm(choice ~ if_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, subset = keep == 1, weight = weight)

texreg::screenreg(list(prepTex(plesm1sr, data = ples), prepTex(plesm1srkeep, data = ples %>% filter(keep == 1)), 
                       prepTex(plesm1sc, data = ples), prepTex(plesm1sckeep, data = ples %>% filter(keep == 1)),
                       prepTex(plesm1fr, data = ples), prepTex(plesm1frkeep, data = ples %>% filter(keep == 1)), 
                       prepTex(plesm1fc, data = ples), prepTex(plesm1fckeep, data = ples %>% filter(keep == 1))))


#1.2.1) Table E1 ####
texreg::texreg(list(prepTex(plesm1srkeep, data = ples %>% filter(keep == 1)), prepTex(plesm1sckeep, data = ples %>% filter(keep == 1)),
                    prepTex(plesm1frkeep, data = ples %>% filter(keep == 1)), prepTex(plesm1fckeep, data = ples %>% filter(keep == 1))),
               file = "./tables/TableE1.tex",
               custom.model.names = c("$I-S$ (Rating)", "$I-S$ (Choice)", "$I-F$ (Rating)", "$I-F$ (Choice)"),
               custom.coef.map = list("is_dist" = "Abs. distance",
                                      "if_dist" = "Abs. distance",
                                      "individual" = "Individual gains",
                                      "societal" = "Societal gains",
                                      "foreign" = "Foreign gains",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Year of Implementation (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = 0.85,
               caption = "Absolute distance and PTA support among respondents with consistent behavior",
               caption.above = T,
               label = "tab:tableE1", 
               custom.note = "\\parbox{0.9\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on respondents. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.11) Interactions announced in the pre-registration but not reported in the paper ####

ples <- ples %>% 
  mutate(lr_cat = factor(ifelse(lrecon < 5, "Left",
                                ifelse(lrecon > 5, "Right",
                                       ifelse(lrecon == 5, "Centre",NA))),
                         levels = c("Left", "Centre", "Right")))

#1.11.1) Figure A1 ####

plesm1sr_lr <- lm(rating ~ is_dist * lr_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sc_lr <- lm(choice ~ is_dist * lr_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fr_lr <- lm(rating ~ if_dist * lr_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fc_lr <- lm(choice ~ if_dist * lr_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)

texreg::screenreg(list(prepTex(plesm1sr_lr, data = ples), prepTex(plesm1sc_lr, data = ples),
                       prepTex(plesm1fr_lr, data = ples), prepTex(plesm1fc_lr, data = ples)))

marginaleffects::avg_slopes(plesm1fr_lr, variable = "if_dist", by = "lr_cat")

predm1srlr <- marginaleffects::predictions(plesm1sr_lr, datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                 lr_cat = c("Left", "Centre", "Right")),
                                           vcov = sandwich::vcovCL(plesm1sr_lr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = lr_cat) %>% 
  mutate(xvar = "IS",
         dv = "R")

predm1sclr <- marginaleffects::predictions(plesm1sc_lr, datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                 lr_cat = c("Left", "Centre", "Right")),
                                           vcov = sandwich::vcovCL(plesm1sc_lr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = lr_cat) %>% 
  mutate(xvar = "IS",
         dv = "C")

predm1frlr <- marginaleffects::predictions(plesm1fr_lr, datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                 lr_cat = c("Left", "Centre", "Right")),
                                           vcov = sandwich::vcovCL(plesm1fr_lr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = lr_cat) %>% 
  mutate(xvar = "IF",
         dv = "R")

predm1fclr <- marginaleffects::predictions(plesm1fc_lr, datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                 lr_cat = c("Left", "Centre", "Right")),
                                           vcov = sandwich::vcovCL(plesm1fc_lr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = lr_cat) %>% 
  mutate(xvar = "IF",
         dv = "C")

rbind(predm1srlr, predm1sclr, predm1frlr, predm1fclr) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper,
                fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line(aes(linetype = mlevel)) + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Ideology") + 
  scale_linetype_discrete("Ideology") +
  theme(legend.position = c(0.2, 0.6))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureA1.pdf",
       width = 20, height = 13.5, units = "cm")

#1.11.2) Figure A2 ####

rbind(predm1srlr, predm1sclr, predm1frlr, predm1fclr) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper,
                fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line(aes(linetype = mlevel)) + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Ideology") + 
  scale_linetype_discrete("Ideology") +
  theme(legend.position = c(0.2, 0.6)) +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureA2.pdf",
       width = 20, height = 13.5, units = "cm")

#1.11.3) Figure A3 ####

plesm1sr_in <- lm(rating ~ is_dist * income_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sc_in <- lm(choice ~ is_dist * income_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fr_in <- lm(rating ~ if_dist * income_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1fc_in <- lm(choice ~ if_dist * income_cat + trade_volume + partner_size + implementation + country, ples, weight = weight)

texreg::screenreg(list(prepTex(plesm1sr_in, data = ples), prepTex(plesm1sc_in, data = ples),
                       prepTex(plesm1fr_in, data = ples), prepTex(plesm1fc_in, data = ples)))

marginaleffects::avg_slopes(plesm1sr_in, variable = "is_dist", by = "income_cat")
marginaleffects::avg_slopes(plesm1fr_in, variable = "if_dist", by = "income_cat")

predm1srin <- marginaleffects::predictions(plesm1sr_in, datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                 income_cat = c("Low", "Med", "High")),
                                           vcov = sandwich::vcovCL(plesm1sr_in, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = income_cat) %>% 
  mutate(xvar = "IS",
         dv = "R")

predm1scin <- marginaleffects::predictions(plesm1sc_in, datagrid(is_dist = seq(from = 0, to = 4, by = .1),
                                                                 income_cat = c("Low", "Med", "High")),
                                           vcov = sandwich::vcovCL(plesm1sc_in, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = is_dist, mlevel = income_cat) %>% 
  mutate(xvar = "IS",
         dv = "C")

predm1frin <- marginaleffects::predictions(plesm1fr_in, datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                 income_cat = c("Low", "Med", "High")),
                                           vcov = sandwich::vcovCL(plesm1fr_in, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = income_cat) %>% 
  mutate(xvar = "IF",
         dv = "R")

predm1fcin <- marginaleffects::predictions(plesm1fc_in, datagrid(if_dist = seq(from = 0, to = 4, by = .1),
                                                                 income_cat = c("Low", "Med", "High")),
                                           vcov = sandwich::vcovCL(plesm1fc_in, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = if_dist, mlevel = income_cat) %>% 
  mutate(xvar = "IF",
         dv = "C")


rbind(predm1srin, predm1scin, predm1frin, predm1fcin) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel, 
                         levels = c("Low", "Med", "High"),
                         labels = c("Low", "Medium", "High")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper,
                fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line(aes(linetype = mlevel)) + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Income") + 
  scale_linetype_discrete("Income") +
  theme(legend.position = c(0.2, 0.6))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureA3.pdf",
       width = 20, height = 13.5, units = "cm")

#1.11.4) Figure A4 ####

rbind(predm1srin, predm1scin, predm1frin, predm1fcin) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("IS", "IF"),
                       labels = c("Individual - Societal", "Individual - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel, 
                         levels = c("Low", "Med", "High"),
                         labels = c("Low", "Medium", "High")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper,
                fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line(aes(linetype = mlevel)) + 
  facet_grid(xvar~., switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Income") + 
  scale_linetype_discrete("Income") +
  theme(legend.position = c(0.2, 0.6))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL)

ggsave(filename = "./figures/FigureA4.pdf",
       width = 20, height = 13.5, units = "cm")

#1.12) Societal gains vs foreign gains ####

#1.12.1) Model 1: Investigation of Absolute Distance ####

plesm1sfr <- lm(rating ~ sf_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sfc <- lm(choice ~ sf_dist + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)

marginaleffects::avg_comparisons(plesm1sfr, variable  = list(sf_dist=c(0, 4)), vcov = sandwich::vcovCL)
marginaleffects::avg_comparisons(plesm1sfc, variable = list(sf_dist=c(0, 4)), vcov = sandwich::vcovCL)

plesm1sfr_cat <- lm(rating ~ sf_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm1sfc_cat <- lm(choice ~ sf_dist_cat + individual + societal +  foreign + trade_volume + partner_size + implementation + country, ples, weight = weight)

#1.12.1.1) Table G1 ####

texreg::texreg(list(prepTex(plesm1sfr, data = ples), prepTex(plesm1sfc, data = ples)),
               file = "./tables/TableG1.tex",
               custom.model.names = c("$S-F$ (Rating)", "$S-F$ (Choice)"),
               custom.coef.map = list("sf_dist" = "Abs. distance",
                                      "individual" = "Individual gains",
                                      "societal" = "Societal gains",
                                      "foreign" = "Foreign gains",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = 0.85,
               caption = "Absolute distance between societal and foreign gains and PTA support",
               caption.above = T,
               label = "tab:tableG1", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on respondents. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.12.2) Absolute Distance and (Dis-)Advantageous Inequality ####

plesm2sfr <- lm(rating ~ sf_dist * sf_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm2sfc <- lm(choice ~ sf_dist * sf_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)

marginaleffects::avg_comparisons(plesm2sfr, variable  = list(sf_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::avg_comparisons(plesm2sfc, variable = list(sf_dist=c(0, 4)), vcov = sandwich::vcovCL, by = "is_dist_bi")

#1.12.2.1) Table G2 ####

texreg::texreg(list(prepTex(plesm2sfr, data = ples), prepTex(plesm2sfc, data = ples)),
               file = "./tables/TableG2.tex",
               custom.model.names = c("$S-F$ (Rating)", "$S-F$ (Choice)"),
               custom.coef.map = list("sf_dist" = "Abs. distance",
                                      "sf_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "sf_dist:sf_dist_biDisadvantageous" = "Abs. dist $\\times$ Disadv. ineq.",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = .85,
               caption = "Absolute distance between societal and foreign gains, (dis-)advantageous inequality and PTA support",
               caption.above = T,
               label = "tab:tableG2", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on individuals. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.12.3) Model 3: Categorical absolute Distance and (Dis-)Advantageous Inequality ####

plesm3sfr <- lm(rating ~ sf_dist_cat * sf_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)
plesm3sfc <- lm(choice ~ sf_dist_cat * sf_dist_bi + trade_volume + partner_size + implementation + country, ples, weight = weight)

marginaleffects::slopes(plesm3sfr, variable = "sf_dist_cat", vcov = sandwich::vcovCL, by = "is_dist_bi")
marginaleffects::slopes(plesm3sfc, variable = "sf_dist_cat", vcov = sandwich::vcovCL, by = "is_dist_bi")

#1.12.3.1) Table G3####

texreg::texreg(list(prepTex(plesm3sfr, data = ples), prepTex(plesm3sfc, data = ples)),
               file = "./tables/TableG3.tex",
               custom.model.names = c("$S-F$ (Rating)", "$S-F$ (Choice)"),
               custom.coef.map = list("sf_dist_catQ2" = "Abs. dist. (Q2)",
                                      "sf_dist_catQ3" = "Abs. dist. (Q3)",
                                      "sf_dist_catQ4" = "Abs. dist. (Q4)",
                                      "sf_dist_biDisadvantageous" = "Disadvantageous inequality",
                                      "sf_dist_catQ2:sf_dist_biDisadvantageous" = "Abs. dist. (Q2) $\\times$ Disadv. ineq.",
                                      "sf_dist_catQ3:sf_dist_biDisadvantageous" = "Abs. dist. (Q3) $\\times$ Disadv. ineq.",
                                      "sf_dist_catQ4:sf_dist_biDisadvantageous" = "Abs. dist. (Q4) $\\times$ Disadv. ineq.",
                                      "trade_volumelarge" = "Trade volume (Large)",
                                      "partner_sizelarge economy" = "Size of partner (Large)",
                                      "implementation2027" = "Implementation Year (2027)",
                                      "countryPL" = "CFE: Poland",
                                      "(Intercept)" = NA),
               float = "htb",
               scalebox = .85,
               caption = "Categorized absolute distance societal and foreign gains, (dis-)advantageous inequality and PTA support",
               caption.above = T,
               label = "tab:tableG3", 
               custom.note = "\\parbox{\\linewidth}{%stars. Entries are unstandardized coefficients from a linear regression model. Standard errors in parentheses are clustered on individuals. Rating is captured on a seven-point scale. Choice is a dummy and we assume linear probabilities.}",
               use.packages = F)

#1.12.4) Preparation Figure G1 ####

predm1sfr <- marginaleffects::predictions(plesm1sfr, newdata = datagrid(sf_dist = seq(from = 0, to = 4, by = .1)),
                                          vcov = sandwich::vcovCL(plesm1sfr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist) %>% 
  mutate(xvar = "SF",
         dv = "R")

predm1sfc <- marginaleffects::predictions(plesm1sfc, newdata = datagrid(sf_dist = seq(from = 0, to = 4, by = .1)),
                                          vcov = sandwich::vcovCL(plesm1sfc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist) %>% 
  mutate(xvar = "SF",
         dv = "C")

predm1sfr_cat <- marginaleffects::predictions(plesm1sfr_cat, newdata = datagrid(sf_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                              vcov = sandwich::vcovCL(plesm1sfr_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist_cat) %>% 
  mutate(xvar = "SF",
         dv = "R",
         mod = "1")

predm1sfc_cat <- marginaleffects::predictions(plesm1sfc_cat, newdata = datagrid(sf_dist_cat = c("Q1", "Q2", "Q3", "Q4")),
                                              vcov = sandwich::vcovCL(plesm1sfc_cat, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist_cat) %>% 
  mutate(xvar = "SF",
         dv = "C",
         mod = "1")

axis_pos <- rbind(ples %>% 
                    group_by(sf_dist_cat) %>% 
                    summarise(pos = mean(sf_dist)) %>% 
                    filter(!is.na(sf_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = sf_dist_cat) %>% 
                    mutate(xvar = "Societal - Foreign")) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("Societal - Foreign")))

bin_1 <- rbind(predm1sfr_cat, predm1sfc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1") %>% 
  filter(dv == "Rating") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar")) %>% 
  mutate(pos = pos * 100)


rbind(predm1sfr, predm1sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  geom_pointrange(data = bin_1, aes(x = pos, y = estimate, ymin = lower, ymax=upper), size = 0.1) +
  facet_grid(xvar~dv, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating")  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p1

bin_1c <- rbind(predm1sfr_cat, predm1sfc_cat) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1") %>% 
  filter(dv == "Choice") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar")) %>% 
  mutate(pos = pos * 100)

rbind(predm1sfr, predm1sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_line() + 
  geom_pointrange(data = bin_1c, aes(x = pos, y = estimate, ymin = lower, ymax=upper), size = 0.1) +
  facet_grid(xvar~dv, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support")  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p1c

#Plot 2
predm2sfr <- marginaleffects::predictions(plesm2sfr, newdata = datagrid(sf_dist = seq(from = 0, to = 4, by = .1),
                                                                        sf_dist_bi = c("Disadvantageous", "Advantageous")),
                                          vcov = sandwich::vcovCL(plesm2sfr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist, mlevel = sf_dist_bi) %>% 
  mutate(xvar = "SF",
         dv = "R",
         mod = "2")

predm2sfc <- marginaleffects::predictions(plesm2sfc, newdata = datagrid(sf_dist = seq(from = 0, to = 4, by = .1),
                                                                        sf_dist_bi = c("Disadvantageous", "Advantageous")),
                                          vcov = sandwich::vcovCL(plesm2sfc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist, mlevel = sf_dist_bi) %>% 
  mutate(xvar = "SF",
         dv = "C",
         mod = "2")

#bins for Plot 2#


predm3sfr <- marginaleffects::predictions(plesm3sfr, newdata = datagrid(sf_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                        sf_dist_bi = c("Disadvantageous", "Advantageous")),
                                          vcov = sandwich::vcovCL(plesm3sfr, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist_cat, mlevel = sf_dist_bi) %>% 
  mutate(xvar = "SF",
         dv = "R",
         mod = "3")

predm3sfc <- marginaleffects::predictions(plesm3sfc, newdata = datagrid(sf_dist_cat = c("Q1", "Q2", "Q3", "Q4"),
                                                                        sf_dist_bi = c("Disadvantageous", "Advantageous")),
                                          vcov = sandwich::vcovCL(plesm3sfc, cluster = ~id)) %>% 
  select(estimate, std.error, lower = conf.low, upper = conf.high, xlevel = sf_dist_cat, mlevel = sf_dist_bi) %>% 
  mutate(xvar = "SF",
         dv = "C",
         mod = "3")

axis_pos <- rbind(ples %>% 
                    group_by(sf_dist_cat) %>% 
                    summarise(pos = mean(sf_dist)) %>% 
                    filter(!is.na(sf_dist_cat)) %>% 
                    arrange(pos) %>% 
                    rename("xlevel" = sf_dist_cat) %>% 
                    mutate(xvar = "Societal - Foreign")) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("Societal - Foreign")),
         pos = pos * 100)

bin_3 <- rbind(predm3sfr, predm3sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 3") %>% 
  filter(dv == "Rating") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar"))

rbind(predm2sfr, predm2sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 1",
         xlevel = xlevel * 100) %>% 
  filter(dv == "Rating") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_pointrange(data = bin_3, aes(x = pos, y = estimate, ymin = lower, ymax=upper, colour = mlevel), size = 0.1, position = position_dodge(width = 10)) +
  geom_line() + 
  facet_grid(xvar~dv, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Rating") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = c(0.3, 0.3),
        legend.text = element_text(size = 8),
        legend.title = element_text(size = 8))  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p2

bin_3c <- rbind(predm3sfr, predm3sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 3") %>% 
  filter(dv == "Choice") %>% 
  dplyr::left_join(., axis_pos, by = c("xlevel", "xvar"))

rbind(predm2sfr, predm2sfc) %>% 
  mutate(xvar = factor(xvar,
                       levels = c("SF"),
                       labels = c("Societal - Foreign")),
         dv = factor(dv,
                     levels = c("R", "C"),
                     labels = c("Rating", "Choice")),
         mlevel = factor(mlevel,
                         levels = c("Disadvantageous", "Advantageous")),
         mod = "Model 2",
         xlevel = 100 * xlevel) %>% 
  filter(dv == "Choice") %>% 
  ggplot(., aes(x = xlevel, y = estimate, ymin = lower, ymax = upper, fill = mlevel)) + 
  geom_ribbon(alpha = 0.2) + 
  geom_pointrange(data = bin_3c, aes(x = pos, y = estimate, ymin = lower, ymax=upper, colour = mlevel), size = 0.1, position = position_dodge(width = 10)) +
  geom_line() + 
  facet_grid(xvar~dv, switch = "y") +
  theme_minimal() +
  labs(x = "Absolute Distance", y = "Predicted PTA Support") +
  scale_fill_discrete("Inequality") + 
  scale_colour_discrete("Inequality") +
  theme(legend.position = "none")  +
  scale_x_continuous(minor_breaks = NULL) +
  scale_y_continuous(minor_breaks = NULL) -> p2c

ggpubr::ggarrange(p1, p2, p1c, p2c)

ggsave(filename = "./figures/FigureG1.pdf",
       width = 20, height = 9, units = "cm")

#1.13) Descriptives ####

#1.13.1) Table C1 #####

quotas_es <- readxl::read_xlsx("./data/quotas_es_pl_gb.xlsx", sheet = 1) 

es_ri <- readRDS("./data/respondent.rds") %>% 
  filter(country == "ES")

es_ri$nuts1 <- factor(ifelse(es_ri$region %in% c("Principado de Asturias", "Cantabria", "Galicia"),
                             "Noroeste",
                             ifelse(es_ri$region %in% c("Arag\u00f3n", "Pa\u00eds Vasco",
                                                        "La Rioja", "Comunidad Foral de Navarra"),
                                    "Noreste",
                                    ifelse(es_ri$region %in% c("Comunidad de Madrid"),
                                           "Com De Madrid",
                                           ifelse(es_ri$region %in% c("Castilla y Le\u00f3n", "Castilla-La Mancha",
                                                                      "Extremadura"),
                                                  "Centro",
                                                  ifelse(es_ri$region %in% c("Islas Baleares", "Catalu\u00f1a",
                                                                             "Comunidad Valenciana"),
                                                         "Este",
                                                         ifelse(es_ri$region %in% c("Andaluc\u00eda", "Regi\u00f3n de Murcia",
                                                                                    "Ciudad Aut\u00f3noma de Ceuta", "Ciudad Aut\u00f3noma de Melilla"),
                                                                "Sur",
                                                                ifelse(es_ri$region %in% c("Canarias"),
                                                                       "Canarias",NA))))))),
                      levels = c("Noroeste", "Noreste", "Com De Madrid", "Centro", "Este", "Sur", "Canarias"))
table(es_ri$nuts1)

es_ri$age_grp <- ifelse(es_ri$age < 35, "< 35",
                        ifelse(es_ri$age < 55, "35 - 54",
                               ifelse(is.na(es_ri$age), NA, "> 54")))

df_es <- bind_rows(es_ri %>% 
                     group_by(nuts1) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = nuts1),
                   es_ri %>% 
                     group_by(gender) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = gender),
                   es_ri %>% 
                     group_by(income_cat) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = income_cat),
                   es_ri %>% 
                     group_by(age_grp) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = age_grp))


dplyr::left_join(quotas_es, df_es) %>% 
  mutate(share_quota = Count / 30,
         share_sample = count_sample / 30) %>% 
  select(Dimension = Dim, Attribute = Attr,
         Quota = Count, "% Quota" = share_quota,
         Sample = count_sample, "% Sample" = share_sample) %>% 
  xtable::xtable(.,
                 digits = c(0, 0, 0, 0, 2, 0, 2),
                 align = c("c", "l", "l", "c", "c", "c", "c"),
                 caption = "Representativeness of the Spanish sample") %>% 
  print(., file = "./tables/TableC1.tex",
        label = "tab:TableC1",
        include.rownames = F, type = "latex", floating	= T, 
        table.placement	= "htb",
        caption.placement = "top")

#1.13.2) Table C2 ####

quotas_pl <- readxl::read_xlsx("./data/quotas_es_pl_gb.xlsx", sheet = 2) 

pl_ri <- readRDS("./data/respondent.rds") %>% 
  filter(country == "PL")

pl_ri$region <- fct_drop(pl_ri$region, only = c("Principado de Asturias", "Cantabria", "Galicia","Arag\u00f3n", "Pa\u00eds Vasco","La Rioja", "Comunidad Foral de Navarra",
                                                "Comunidad de Madrid","Castilla y Le\u00f3n", "Castilla-La Mancha","Extremadura",
                                                "Islas Baleares", "Catalu\u00f1a", "Comunidad Valenciana","Andaluc\u00eda", "Regi\u00f3n de Murcia","Ciudad Aut\u00f3noma de Ceuta", 
                                                "Ciudad Aut\u00f3noma de Melilla","Canarias"))

levels(pl_ri$region) -> regions


pl_ri$nuts1 <- factor(ifelse(pl_ri$region %in% regions[c(3,13)],
                             "Centralny",
                             ifelse(pl_ri$region %in% regions[c(6,12)],
                                    "Poludniowy", 
                                    ifelse(pl_ri$region %in% regions[c(4,9,10)],
                                           "Wschodni",
                                           ifelse(pl_ri$region %in% regions[c(16,17,5)],
                                                  "Polnocno-Zachodni", 
                                                  ifelse(pl_ri$region %in% regions[c(1,8)],
                                                         "Poludniowo-Zachodni", 
                                                         ifelse(pl_ri$region %in% regions[c(2,14,11)],
                                                                "Polnocny", 
                                                                ifelse(pl_ri$region %in% regions[c(15,7)],
                                                                       "Wojewodztwo Mazowieckie",NA))))))),
                      levels = c("Centralny", "Poludniowy", "Wschodni", "Polnocno-Zachodni", "Poludniowo-Zachodni", "Polnocny", "Wojewodztwo Mazowieckie"))


table(pl_ri$nuts1)

pl_ri$age_grp <- ifelse(pl_ri$age < 35, "< 35",
                        ifelse(pl_ri$age < 55, "35 - 54",
                               ifelse(is.na(pl_ri$age), NA, "> 54")))

df_pl <- bind_rows(pl_ri %>% 
                     group_by(nuts1) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = nuts1),
                   pl_ri %>% 
                     group_by(gender) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = gender),
                   pl_ri %>% 
                     group_by(income_cat) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = income_cat),
                   pl_ri %>% 
                     group_by(age_grp) %>% 
                     summarise(count_sample = n()) %>% 
                     rename(Attr = age_grp))

dplyr::left_join(quotas_pl, df_pl) %>% 
  mutate(share_quota = Count / 30,
         share_sample = count_sample / 30) %>% 
  select(Dimension = Dim, Attribute = Attr,
         Quota = Count, "% Quota" = share_quota,
         Sample = count_sample, "% Sample" = share_sample) %>% 
  xtable::xtable(.,
                 digits = c(0, 0, 0, 0, 2, 0, 2),
                 align = c("c", "l", "l", "c", "c", "c", "c"),
                 caption = "Representativeness of the Polish sample") %>% 
  print(., file = "./tables/TableC2.tex",
        label = "tab:TableC2",
        include.rownames = F, type = "latex", floating	= T, 
        table.placement	= "htb",
        caption.placement = "top")

#1.13.1) Table C3 #####

bind_rows(ples %>% 
            select(country, rating, choice, is_dist, if_dist, sf_dist, individual, societal, foreign,
                   weight, income, lrecon) %>% 
            filter(country == "ES") %>% 
            psych::describe(.) %>% 
            data.frame() %>% 
            mutate(country = "ES"),
          ples %>% 
            select(country, rating, choice, is_dist, if_dist, sf_dist, individual, societal, foreign,
                   weight, income, lrecon) %>% 
            filter(country == "PL") %>% 
            psych::describe(.) %>%
            data.frame() %>% 
            mutate(country = "PL")) %>% 
  mutate(vars = factor(vars,
                       levels = c(1:12),
                       labels = c("Country", "Rating", "Choice", "Abs. distance (I-S)", "Abs. distance (I-F)", "Abs. distance (S-F)",
                                  "Individual", "Societal", "Foreign",
                                  "Weight", "Income", "Ideology"))) %>% 
  filter(vars != "Country") %>% 
  select(Country = country, Variable = vars, n, Mean = mean, SD = sd, Min = min, Max = max) %>% 
  xtable::xtable(.,
                 digits = c(0,0, 0, 0, 2, 2, 2, 2),
                 align = c("c", "c", "l", "c", "c", "c", "c", "c"),
                 caption = "Descriptive statistics") %>% 
  print(., file = "./tables/TableC3.tex",
        label = "tab:TableC3",
        include.rownames = F, type = "latex", floating	= T, 
        table.placement	= "htb",
        caption.placement = "top")
